#include <navierstokes.hpp>


#define WEAKCL 1

namespace Life
{


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::updatePtFixe(element_fluid_type & U, element_fluid_type & Uold,element_fluid_velocity_type & mesh_velocity,ale_map_element_type & depl )
    {

        using namespace Life::vf;

        //space_fluid_ptrtype XhBis = U.functionSpace();

        vector_ptrtype F( M_backend->newVector( Xh ) );
        sparse_matrix_ptrtype A( M_backend->newMatrix(Xh,Xh) );

        mesh_ptrtype mesh = Xh->mesh();

        element_fluid_type V( Xh, "v" );
        element_fluid_0_type u = U.template element<0>();
        element_fluid_0_type v = V.template element<0>();
        element_fluid_1_type p = U.template element<1>();
        element_fluid_1_type q = V.template element<1>();

        element_fluid_0_type uold = Uold.template element<0>();
        element_fluid_1_type pold = Uold.template element<1>();

        //--------------------------------------------------------------------------------------------------//

        double EE = 21*1e5;
        double sigma = 0.28;
        double mu = EE/(2*(1+sigma));
        double lambda = EE*sigma/((1+sigma)*(1-2*sigma));

        AUTO( Iddd , oneX()*trans(oneX()) + oneY()*trans(oneY()) );
        AUTO( Fv , Iddd + gradv(depl) );
        AUTO( Ev, 0.5*(gradv(depl)+trans(gradv(depl)) ) + 0.5*trans(gradv(depl))*gradv(depl) );
        AUTO( Sv, lambda*trace(Ev)*Iddd + 2*mu*Ev );

        //--------------------------------------------------------------------------------------------------//

        //Tenseur des deformations
        AUTO( deft, 0.5*(gradt(u)+trans(gradt(u)) ) );
        AUTO( def, 0.5*(grad(v)+trans(grad(v)) ) );
        //Matrice Identitee
        AUTO (Id , oneX()*trans(oneX()) + oneY()*trans(oneY()) );
        // Tenseur des contraintes (trial)
        AUTO(Sigmat, -idt(p)*Id + 2*mu*deft );
        // Tenseur des contraintes (test)
        AUTO(Sigma, -id(q)*Id + 2*mu*def );
        //force volumique :
        AUTO(f,-rho*cst(9.81)*vec(cst(0.),cst(1.)) );

        //AUTO(vitDepl, vec(cst(1.),cst(0.) ) );
        //AUTO(vitDepl, vec(cst(0.),-0.3*sin(M_PI*Px()/4.) ) );
        //AUTO(vitDepl, -0.3*sin(M_PI*Px()/4.)*N() );

        // Vitesse d'entree du fluide
        AUTO(Vit_Entree, 10*Vitesse_Entree*( Py()*(1.0-Py()) )*oneX() );


        //--------------------------------------------------------------------------------------------------//
        //terme non lineaire
        form2( Xh,Xh, A, _init=true )  =
            integrate (elements(mesh), _Q<3*Order>(),
                       rho*trans(id(v))*(gradt(u))*(idv(uold) - idv(mesh_velocity)) );

        // sigma : grad(v) sur Omega
        form2( Xh, Xh, A ) +=
            integrate( elements(mesh), _Q<2*(Order)>(),
                       + trace(Sigmat*trans(grad(v))) );

        // terme assurant l'incompressibilite
        form2( Xh, Xh, A ) +=
            integrate( elements(mesh), _Q<2*(Order-1)>(),
                       + divt(u)*id(q) );

        // terme necessaire pour condition forte avec petsc
        form2( Xh, Xh, A ) +=
            integrate( elements( mesh ),// _Q<3*uOrder-1+3*(GeoOrder-1)>(),
                       +0*div(v)*idt(p)+0*divt(u)*id(q) +0*idt(p)*id(q) );

        /*
        //terme en plus pour div=0 (voir these goncalo)
        form2( Xh, Xh, A ) +=
        integrate( elements( mesh ),// _Q<3*uOrder-1+3*(GeoOrder-1)>(),
        divv(uold)*trans(idt(u))*id(v) );
        */

        // forces volumique
        form1( Xh, F, _init=true ) =
            integrate( elements(mesh), _Q<3*(Order)>(),
                       trans(f)*id(v) );

        //imposition de la contrainte
        form1( Xh, F ) +=
            integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<3*Order>(),
                       + trans(Fv*Sv*N())*id(v) );
        form1( Xh, F ) +=
            integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<3*Order>(),
                       + trans(Fv*Sv*N())*id(v) );

        //--------------------------------------------------------------------------------------------------//

#if WEAKCL
        // terme de bord classique
        form2( Xh, Xh, A ) +=
            integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                       - trans(Sigmat*N())*id(v) );
        //+ trans(idv(uold)*trans(idt(u))*N())*id(v) );
        form2( Xh, Xh, A ) +=
            integrate( markedfaces(mesh,mesh->markerName("Paroi")), _Q<2*Order>(),
                       - trans(Sigmat*N())*id(v) );
        //+ trans(idv(uold)*trans(idt(u))*N())*id(v) );

        // pas de terme car contrainte imposee
        /*
          form2( Xh, Xh, A ) +=
          integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
          - trans(Sigmat*N())*id(v) );
          //+ trans(idv(uold)*trans(idt(u))*N())*id(v) );
          form2( Xh, Xh, A ) +=
          integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
          - trans(Sigmat*N())*id(v) );
          //+ trans(idv(uold)*trans(idt(u))*N())*id(v) );
          */

        //terme de bord weak boundary
        form2( Xh, Xh, A ) +=
            integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                       - trans(Sigma*N())*idt(u) // symetrisation
                       //+ trans(idv(uold)*trans(id(v))*N())*idt(u)
                       + penalbc*trans(idt(u))*id(v)/hFace() );
        form2( Xh, Xh, A ) +=
            integrate( markedfaces(mesh,mesh->markerName("Paroi")), _Q<2*Order>(),
                       - trans(Sigma*N())*idt(u) //symetrisation
                       //+ trans(idv(uold)*trans(id(v))*N())*idt(u)
                       + penalbc*trans(idt(u))*id(v)/hFace() );
        form2( Xh, Xh, A ) +=
            integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
                       //- trans(Sigma*N())*idt(u) //pas de symetrisation si contrainte imposee!!
                       //+ trans(idv(uold)*trans(id(v))*N())*idt(u)
                       + penalbc*trans(idt(u))*id(v)/hFace() );
        form2( Xh, Xh, A ) +=
            integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
                       //- trans(Sigma*N())*idt(u) //pas de symetrisation si contrainte imposee!!
                       //+ trans(idv(uold)*trans(id(v))*N())*idt(u)
                       + penalbc*trans(idt(u))*id(v)/hFace() );


        form1( Xh, F ) +=
            integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                       - trans(Sigma*N())*Vit_Entree // symetrisation
                       //+ trans(idv(uold)*trans(id(v))*N())*Vit_Entree
                       + penalbc*trans(Vit_Entree)*id(v)/hFace() );

        form1( Xh, F ) +=
            integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
                       //  - trans(Sigma*N())*idv(mesh_velocity) //pas de symetrisation si contrainte imposee!!
                       //+ trans(idv(uold)*trans(id(v))*N())*idv(mesh_velocity)
                       + penalbc*trans(idv(mesh_velocity))*id(v)/hFace() );
        form1( Xh, F ) +=
            integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
                       //- trans(Sigma*N())*idv(mesh_velocity) //pas de symetrisation si contrainte imposee!!
                       //+ trans(idv(uold)*trans(id(v))*N())*idv(mesh_velocity)
                       +penalbc*trans(idv(mesh_velocity))*id(v)/hFace() );

#endif


        //termes temporels (ordre 1 pour le moment)
        form2( Xh, Xh, A ) +=
            integrate( elements(mesh), _Q< 2*Order >(),
                       + trans(idt(u))*id(v)*M_bdf_fluid->polyDerivCoefficient(0));
        element_fluid_type Buzz( Xh, "v" );
        element_fluid_0_type buzz = U.template element<0>();
        buzz.container() = M_bdf_fluid->polyDeriv().template element<0>().container();
        form1( Xh, F ) +=
            integrate( elements(mesh), _Q< 2*Order >(),
                       + trans(idv(buzz))*id(v) );

#if 1

        //AUTO( beta_N , trans(N())*idv(beta) );
        AUTO( beta_abs , vf::sqrt( trans(idv(uold)-idv(*M_mesh_velocity))*(idv(uold)-idv(*M_mesh_velocity)) ));

        //terme de stabilisation
        form2( Xh, Xh, A ) +=
            integrate( internalfaces(mesh), _Q< 3*Order >(),
                       beta_abs*vf::pow(hFace(),2.0)*(1./vf::pow(cst(double(Order)),3.5))*trans(jumpt(gradt(u)))*jump(grad(v)) );

#endif

        A->close();
        F->close();

#if !WEAKCL
        form2( Xh, Xh, A ) +=
            on( markedfaces(mesh, mesh->markerName("Paroi")) , u, F, vec(cst(0.),cst(0.)) );
        form2( Xh, Xh, A ) +=
            on( markedfaces(mesh, mesh->markerName("ParoiH")) , u, F, idv(mesh_velocity) );
        form2( Xh, Xh, A ) +=
            on( markedfaces(mesh, mesh->markerName("ParoiB")) , u, F, idv(mesh_velocity) );
        form2( Xh, Xh, A ) +=
            on( markedfaces(mesh, mesh->markerName("Entree")) , u, F, Vit_Entree );
#endif

        vector_ptrtype X( M_backend->newVector( Xh ) );
        M_backend->solve( _matrix=A, _solution=X, _rhs=F );
        U = *X;

    }


} //end namespace Life


template class Life::NavierStokes<2,2,Life::Simplex>;
